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Abstract 

Model of laminated wave turbulence allows to study statistical 
and discrete layers of turbulence in the frame of the same model. 
Statistical layer is described by Zakharov-Kolmogorov energy spectra 
in the case of irrational enough dispersion function. Discrete layer 
is covered by some system (s) of Diophantine equations while their 
form is determined by wave dispersion function. This presents a very 
special computational challenge - to solve Diophantine equations in 
many variables, usually 6 to 8, in high degrees, say 16, in integers of 
order 10 16 and more. Generic algorithms for solving this problem in 
the case of irrational dispersion function have been presented in our 
previous papers. In this paper we present a new algorithm for the case 
of rational dispersion functions. Special importance of this case is due 
to the fact that in wave systems with rational dispersion the statistical 
layer does not exist and the general energy transport is governed by 
the discrete layer alone. 
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1 Introduction 



The general theory of fluid mechanics begins in 1741 with the work of Leon- 
hard Euler who was invited by Frederick the Great to construct an intrinsic 
system of water fountains. Euler began with deducing the equations which 
are now called Euler equations; they describe the ideal (inviscid) liquid and 
are derived from the classical Newton's conservation laws written for a fluid 
particle. Euler equations, regarded with various boundary conditions and 
specific values of some parameters describe enormous number of wave sys- 
tems, for instance, capillary waves, surface water waves, atmospheric plan- 
etary waves, drift waves in plasma, Tsunami, freak waves, etc. The general 
form of reduced Euler equations suitable for studying one specific type of 
waves can be written as 



where £ and TV are linear and nonlinear operators correspondingly, and e is a 
small parameter chosen according to the properties of the wave system under 
consideration. For instance, it can be taken as a ratio of wave amplitude to 
its length, or as a ratio of a particle velocity to the phase velocity, or some 
other way. A linear wave is then a solution of the corresponding linear 
equation C{ip) = and has standard form Aexpi[kx — uit] with amplitude 
A, wave vector k and dispersion function Ui = uj(ki). The form of dispersion 
function is defined then by boundary conditions. The existence of a small 
parameter e allows to reduce the study of all nonlinear waves to those which 
are resonantly interacting, that is, satisfy resonant conditions 



Notice that amplitudes of resonantly interacting waves are not constant any 
more and standard multi-scale method yields the corresponding system of 
ordinary differential equations (ODEs) on these amplitudes. The energetic 
behavior of a wave system depends drastically on whether wave vectors fcj 
have real or integer coordinates. The first case (real-valued coordinates) is 
treated in the frame of statistical wave turbulence (SWT) theory pQ, with 
additional assumption that u(ki)/u(kj) is an algebraic number of degree > 2. 
The energy transport in these systems is covered by the wave kinetic equa- 
tion. The second case (integer-valued coordinates) is described by discrete 
wave turbulence (DWT) theory [2J, and energy transport is presented by 
a few quasi-periodic processes. Model of laminated turbulence |3] presents 
SWT and DWT as two layers of a wave system, with elaborate transition 
from one layer to another. One of the novel problems emerging from this 
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model is the necessity to solve ([T]) for very big integers. 



In the first two articles [I], [5] of this series we presented algorithms for 
finding resonant wave interactions for irrational dispersion functions, with 
two illustrative examples: (1) gravitational water waves, uo = \Jm 2 + n 2 (4- 
wave interactions); and (2) ocean planetary waves, u = 1/ ' \/m 2 + n 2 (3- wave 
interactions). The key points of the presentation were, first, that our algo- 
rithms for these cases differ only in some details and their core is applicable 
to a wide class of dispersion functions, thus justifying the name of "generic". 
Second, irrational equations in integers were solved without use of floating- 
point arithmetic and not even resolving the irrationalities involved. This 
gave us an enormous gain both in performance time and orders of numbers 
used. 

In the present paper we construct a special algorithm for solving (CQ) in 
case of a rational dispersion function. Notice that for any rational dispersion 
function, u(ki)/u(kj) is obviously a rational number, that is, an algebraic 
number of degree 1. It makes SWT theory not applicative for these type of 
wave systems because statistical layer of turbulence does not exists and the 
whole energetic behavior is covered by the discrete layer only. This makes 
the creation of some fast algorithm for computing integer solutions of ([I]) for 
the case of rational dispersion function of high importance. 

2 General idea of the algorithm 

Obviously, any equation in rational functions in integers can be trivially 
transformed into a Diophantine equation. For 



which, however, leads to huge powers and extensive search. The idea un- 
derlying our algorithm is quite simple and we illustrate it by the example 
below. 




(2) 



the corresponding Diophantine equation will be 
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Example Suppose we need solve in integers an equation 

P 

a = b—, < a < do, < b < 6 



(4) 



where P/Q is an irreducible fraction. We could transform it into aQ = bP 
and perform exhaustive search in the region 0<a<ao,0<6<6o with 
computational complexity O(a bo). 

However, we notice that the number b^ is integer only if b is a multiple 
of the denominator Q. Then (a, b) is a solution only if b = kQ with integer 
k. Which immediately gives a = kP and (kP, kQ) is a solution for any 
k, 1 < k < min(P/ao, Q/bo) and these are all the solutions of the equation. 
Notice that there is no search at all involved. 

To show the power of the approach outlined above in practice, we proceed 
further with the example of spherical planetary waves. 

2.1 Example 1: spherical planetary waves 

The turbulence of the spherical planetary waves is governed by the barotropic 
vorticity equation on a sphere 



dt 



where 



d 1 ^ 
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cos 2 d\ 2 
A linear wave has the form 



tan. 
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(5) 



and J (a, b) 



COS( 



,da db 
dXd^p 



da db , 
dSdX' 



AP™ (sin 4>) exp i[mX + 



2m 



n(n + 1) 



it"- 



with constant wave amplitude A, dispersion function uj = —2m/[n(n + 1)] 
and P™(x) being the associated Legendre function of degree n and order 
m. Resonance conditions in this case have form[6]: 



m l + m 2 — m 3 

rrii < rti Vi = 1, 2, 3 
\n% -n 2 \<nj,<ni+n2 
n\ + 112 + n-s = 1( mod 2) 



(6) 
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where u) { = 7714/(7^(7^+1)) . 

We are going to find all the solutions of Sys. (E|) in a finite domain D, 
i.e. < mi,ni < D Vz = 1,2,3. In our numerical experiments we operated 
with D = 1000 further called the main domain. 

2.2 Computational Preliminaries 

The straightforward approach would be to multiply the first equation of Sys. 
(jEJ) with all three denominators rij (77^+1), substitute 777,3 with mi + m 2 and 
perform full search on 777,1, m 2 , 77,1, 77,2, 713. This evidently implies .D 5 compu- 
tation time and operating with numbers of the order of D 5 . For the main 
domain D = 1000 this is halfway feasible with a large computer but clearly 
not for everyday use with a usual PC. Moreover, computation time and order 
of numbers used grow rapidly with the domain, so when need for computa- 
tions in larger domains arises, as it surely will, the algorithm will fail. 

We are going to present a far more efficient algorithm. 

2.3 Algorithm Description 

• Step 1: Search on rii, n 2 , n 3 . 

The search on 77,1,77,2,77,3 is organized conventionally. Without loss of 
generality, consider 711 < 77,2. Notice that 77,3 always lies between n\ and 
rz.2 and ^from the two "triangle inequalities" of Sys. ([6]) the second one 
always holds, while the first one implies n 3 > n 2 — rti which limits the 
search on n 3 if 77-2 — 72! > n 1 . The oddity condition allows us to run the 
cycle on n 3 in steps of 2. 

Up to now, the computational complexity is 0(D 3 ). 

• Step 2: Cycles elimination on mi,m 2 . 

The numbers of the form n(n + 1) are sometimes called "box numbers" 
(analogous to the square numbers n 2 ) and we introduce notation 6, = 
riiijii + 1). Now we rewrite the first equation of Sys. ([6]) as 

7711/&1 + m 2 /b 2 = rrii/b 3 + m 2 /b 3 (7) 
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or 

^3 ~ h b 2 

m 2 = m 1 - — • — (8) 

b 2 - h h 

Let us find the greatest common divisor GCD of the numerator and 
denominator of the fraction on the right side and reduce by it. The 
equation now has the form 

m 2 =m 1 - — (9) 
tin 

and every solution has the form 

mi = kR D , m 2 = kR N , k < min{rii/ R D) n 3 / (R N + R D )). (10) 
The second condition follows from 

m 3 = mi + m 2 < n 3 (11) 
and is stronger than m 2 < n 2 . 

The computational complexity of the whole algorithm is thus 0(log-D-D 3 ) 
D 3 for the cycle on rii and logD for the GCD. 



Remark The algorithm above implies operating with numbers of the 
order of D 4 in one certain place, namely, transforming 

h ~ h h Rn , x 

b 2 -h bx R D 1 ] 

This could lead to overflows be D large and computer small, say D = 
1000 and 32 bit computer or D = 10 6 and 64 bit computer. There 
is, however, an elegant way to avoid difficulties at this point which we 
describe in the next Step. 



• Step 3: Avoiding multiplications. 



Given the fraction product, we first reduce 63 — b\ and b 2 — 63 by their 
GCD, then b 2 and b\ by their GCD. This leaves us with a product of 
two irreducible fractions (r 3 i/r 23 ) ■ (r 2 /r 1 ). Now we reduce crosswise: 
r 3 i and r 1; r 23 and r 2 . The last reduction gives an "irreducible product" 
of two fractions {rr 3 i/rr 23 ) ■ {rr 2 /rr\), i.e. had we performed the multi- 
plications, the resulting fraction would stay irreducible. The reduction 
schema is presented in Fig. [U 
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Figure 1: Bringing a product of two fractions to complete irreducibility with- 
out multiplying. 

We still do not perform multiplications for fear of an overflow. But 
now it is evident that a solution can only exist if rr 23 < rii, rr\ < 
n>i, rr 31 < n,2, rr2 < ri2- We first check these inequalities; if one or 
more of them do not hold, we proceed with the n-cycle, otherwise we 
may safely perform multiplications (both products do not exceed D 2 ) 
and look for solutions. 

2.4 Example 2: drift waves in a channel 

The turbulence of the drift waves is described by the same equation as in 
Sec. 2.1 but in Descartes coordinates and in the infinite channel [7J. In this 
case dispersion function has a (slightly simplified) form u = 2m/ (n 2 + 1) and 
resonance conditions are 

' Ui + to 2 = UJ3 
^ m 1 +m 2 =m 3 
rrii <rii Vi = 1, 2, 3 
Ui^rij Vi^j 

Search cycles on rii, n 2 , n 3 become somewhat more extensive due to the 
lack of the two conditions mentioned above. On the other hand, the core of 
the algorithm - the four reductions (Step 2) - are preserved one-to-one, as 
well as the post- reduction overflow check (Step 3). 

The computational complexity of the whole algorithm is also 0(\ogDD 3 ) 
as in the previous case. 
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3 Numerical results and some discussion 



Our algorithm has been implemented in VBA programming language; for 
D = 1000 computation time (without disk output of solutions found) on a 
low-end PC (800 MHz Pentium III, 512 MB RAM) is about 7.5 minutes. 
Altogether 7282 solutions (example 1) have been found. Some overall nu- 
merical data is given in the Tables and Figures below. 
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Figure 2: Example 1: Number of all solutions in partial domains: squares 
(points with diamonds) and circles (points with circles). 

In FigfS] the number of solutions in partial domains is shown for the first 
example (atmospheric planetary waves) and we conclude that the solutions 
are concentrated along X and Y axes. In Figj3] the histogram of vector mul- 
tiplicities is presented which shows in how many solutions one vector can 
participate. On the axis X the multiplicity of a vector is shown and on the 
axis Y the number of vectors with a given multiplicity. One can see immedi- 
ately that most part of vectors take part only in one solution and multiplicity 
decreases exponentially with the number of solutions. 

As to our second example (drift waves in a channel) we notice, first of all, 
much less solutions (477) in the same main domain D = 1000. Therefore, 
not much can be said about the asymptotic of solution number in partial 
domains. There is no need to present multiplicities graphically in this case. 
In the whole calculation domain D = 1000 there is just one vector (1, 5) 
participating in solutions with multiplicity 5, one vector (78, 99) with multi- 
plicity 4 and the overall distribution is as follows: 
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Figure 3: Example 1: Histogram of vector multiplicities. 



Table 1. Example 2: Vector multiplicities. 



In order to understand the energetic behavior of 2-dimensional discrete 
wave system, the standard way is to present is graphically on the integer 
lattice in following way. Each node with coordinates m, n presents a corre- 
sponding wave vector k = (m, n) and nodes- vectors are connected by lines is 
they are parts of the same solution. An example of this geometrical structure 
is given in FigJU 



50 100 150 200 250 300 



Figure 4: Example of geometrical structure of a solution set. 



This geometrical representation is needed in order to understand what 
sort of equations (ODEs) for the amplitudes of resonantly interacting waves 
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we have to solve. Namely, one single triangle in /c-space corresponds to 

Ai = a x A 2 A 3 

A 2 = a 2 A x A 3 (14) 
A 3 = a 3 A 1 A 2 

where coefficients a« are known functions on m^rij, % = 1,2,3. If one wave 
takes part in two solutions, we get two systems of this form connected via 
this wave, for instance, if the second solution corresponds to 

A A = a 4 A 5 A 6 , A 5 = a 5 A 4 A 6 , A 6 = a 6 A 4 A 4 , 

and they are connected via one wave, say A 3 = At, then corresponding 
system of ODEs takes form 



'A, 


= a x A 2 A 3 


A 2 


= a 2 A x A 3 


< A 3 


= \{a 3 AxA 2 + a 4 A 5 A 6 ) 


is 


= a 5 A 3 A 6 


A 


= a 6 A 3 A 5 



and so on. 

Obviously, the geometrical structure is too confusing to be informative 
and what we really need is a topological structure of a solution set, i.e. the 
graph formed by triangles as primary elements. Namely, it is enough to 
compute all non-isomorphic topological elements because all isomorphic ele- 
ments are described by the same system of ODEs. For instance, all primary 
elements (isolated resonant triads) are described by (|14p . all "butterflies" 
(groups of two connected triads) are described by (TT5|) . etc. The difference 
between two isomorphic topological elements lies in the coefficients a» which 
are functions of the wave numbers, and therefore, will take different magni- 
tudes for different resonant triads. Topological structure of the solution set 
for our first example is shown in Fig. 0for domain m,n < 50. This domain 
contains 42 solutions: 15 isolated triangles, two "butterflies" (groups of two 
connected triangles), one chain of three connected triangles and two more 
complex graphs. 



4 Summary 

This paper concludes the series of three papers on generic algorithms for lam- 
inated wave turbulence. We have presented algorithms for polynomial disper- 
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Figure 5: Example 1: Topological structure of solutions for D = 50. 



sion function depending irrationally on the wave vector length k = y/m 2 + n 2 
and for an arbitrary rational dispersion function. We have also shown that 
the topological elements of the solution set (for the discrete layer of turbu- 
lence) give the whole information about the energy transport in these wave 
systems. In fact we applied this approach for our first example (atmospheric 
planetary waves) and studied all the topological elements in the meteorolog- 
ically significant domain (m,n < 21) for climate range processes[8j. More 
precisely, we have found analytically solutions of corresponding systems of 
the form (|14p in terms of Jacobean elliptic functions and computed their pe- 
riods and other properties for characteristic meteorological data. As a result, 
a novel model of the known physical phenomena - intra-seasonal oscillations 
in the Earth atmosphere - has been developed. 

Our further interest lies now in the area of symbolic computations. In- 
deed, beginning with an equation like (jSj) with given boundary conditions, 
we have a completely constructive procedure of obtaining precise form for: I) 
resonance conditions (JTJ); II) coefficients a, of the primary element ( 1141) ; III) 
all topological elements - as graphs and as corresponding systems of ODEs. 
All this can be programmed symbolically in MATHEMATICA and solutions 
can be found using our generic algorithms. This is an on-going work now and 
we plan to create a useful program tool for making basic research in the area 
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of discrete wave turbulence. There are some open mathematical questions 
yet to be solved - for instance, the problem of graph isomorphism appearing 
at the step when all different topological elements have to be computed. 

Another possible development would be the study of the 4-wave interac- 
tions, that is, with primary elements being not triads but quartets of waves 
(see example of gravitational water waves [I], [5] )• The same constructive pro- 
cedure as for 3-wave interactions can be applied but the resulting topology 
will be much more complicated due to the principal difference between 3- 
and 4-wave systems. In 3-wave system there exist the only mechanism for 
the energy transport - transport over the scales. In 4-wave system there are 
two qualitatively different mechanisms of the energy flow - over the scales 
and over the phases, and they can combine in a highly nontrivial way. 

Acknowledgement. E.K. acknowledges the support of the Austrian 
Science Foundation (FWF) under projects SFB F013/F1304. 
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